%% Jicamarca ISR data processing
% The script was written to comapre the drift at 150 and 600 km
% the original file is from P. Alken

%# Field 1: timestamp (UT)
%# Field 2: local time (hours)
%# Field 3: year
%# Field 4: month
%# Field 5: day
%# Field 6: altitude (km)
%# Field 7: v_z (m/s)
%# Field 8: sigma (m/s)
%# Field 9: E (mV/m)
%# Field 10: sigma_E (mV/m)
%# Field 11: B (nT)


load /nfs/satmag_work/mnair/projects/eef/jicamarca_isr.dat
k = ones([1,length(jicamarca_isr)])';
fday_jic = datenum([1970.*k 1.*k 1.*k 0.*k 0.*k jicamarca_isr(:,1)]) ;
ntimestamps = unique(fday_jic);

%%
for i = 1: length(ntimestamps),
    
    L = fday_jic == ntimestamps(i);
    
    driftdata = jicamarca_isr(L,7);
    sigma = jicamarca_isr(L,8);
    
    LL = jicamarca_isr(L,6) >=100 & jicamarca_isr(L,6) <= 300 & sigma < 1;
        
    jicamarca_drift_150(i) = nanmean(driftdata(LL));
    jicamarca_sigma_150(i) = nanmean(sigma(LL));
    jicamarca_std_150(i) = std(driftdata(LL));
    
    LL = jicamarca_isr(L,6) >=500 & jicamarca_isr(L,6) <= 700 & sigma < 1;
    
    jicamarca_drift_600(i) = nanmean(driftdata(LL));
     jicamarca_sigma_600(i) = nanmean(sigma(LL));
    jicamarca_std_600(i) = std(driftdata(LL));
    
   
    
    LL = jicamarca_isr(L,6) >=500 & sigma < 1;
    
    jicamarca_drift_gt_500(i) = nanmean(driftdata(LL));
     jicamarca_sigma_gt_500(i) = nanmean(sigma(LL));
    jicamarca_std_gt_500(i) = std(driftdata(LL));
    
     
    % All Altitudes with sigma < 1
    LL = sigma < 1;
    
    jicamarca_drift(i) = nanmean(driftdata(LL));
    jicamarca_sigma(i) = nanmean(sigma(LL));
    jicamarca_std(i) = std(driftdata(LL));
    
    % All Altitudes
    
    
    jicamarca_drift(i) = nanmean(driftdata);
    jicamarca_sigma(i) = nanmean(sigma);
    jicamarca_std(i) = std(driftdata);
    
   
    % Jicamarca Radar Measurements are most accurate for the
    % interval 200 - 400 km altitude, according to several studies. S
    % otrying it out
    
    LL = jicamarca_isr(L,6) >=200 & jicamarca_isr(L,6) <= 400;
    
    jicamarca_drift_500(i) = nanmean(driftdata(LL));
     jicamarca_sigma_500(i) = nanmean(sigma(LL));
    jicamarca_std_500(i) = std(driftdata(LL));

end

%Results
%=========
% The script was made to comapre the drift parameters at two heights
% drift velocity - very similar during local day time. Differing in night with the
% averages at 600 km more negative than the averages at 200. But when the
% drift is almost near zero in the night (storm ?), both the velocities
% match

% sigma (is this the error ?). 150 less all the time than 600. The sigma at
% 600 is highest around morning 7 am and lowest during local noon. sigma ta
% 150 (or 200 :-) ) doesn't show this variation - flat

% std of the averages. Standard deviation of the data. During  night time,
% the lower drift (200) shows considerabli higher std that that
% of the drift around 600 km. During day time, most;ly they are much
% smaller and similar. 

 save /data/backup/mnair/longp/jicamarca_isr_2000_2005_sigma ;
%%
% Use Jicamarca day-night data sent by Anderson . Create the data_index
% with equal interval for spectral analysis.

load ('/data/backup/mnair/ace/jcamarca_isr_day_night.mat');

min_time_length = 6/24;
lt_start = 8;
lt_end = 16;
tol = 0.004;
np =0;
nd =1;

fday_lt = ((fday - 5/24) - floor (fday - 5/24)) * 24;

for i = 1: 25464,
    
    if (fday(i+1) - fday(i) <= tol & fday_lt(i) > lt_start & fday_lt(i) < lt_end ...
            & fday(i) - fday(i - np) <= min_time_length )
        np = np + 1;
    else
       % if (np >=16)
        if ( fday(i) - fday(i - np ) >= min_time_length )
            data_index(nd,2) = i;
            data_index(nd,1) = i - np;
            np = 0;
            nd = nd + 1;
        else
            np=0;
            
        end;
    end;
    
end;

% The processing was done on the wondows computer
fejer= load('/data/backup/mnair/longp/jic_isr_fejer.out');
fejer(end,:) = [];% For some reason the model duplicated the last point
% now creating a data matrix similar to that of eef supplied by
% alken. Note that some columns with 0s will be there
% mutiplying the drift with 23000 nT for getting eef

eef(:,1) = fday - datenum(2000,1,1);
eef(:,4) = fday_lt;
eef(:,6) = drift*23000/1e6;
eef(:,15) = fejer(:,2)*23000/1e6;

save /data/backup/mnair/ace/jcamarca_isr_fejer.mat
%% Remove the climatological variation

% write the data to a file for running Fejer's model

% load f107 data 
load /data/backup/mnair/indices/f107.dat;
load ('/data/backup/mnair/ace/jcamarca_isr_day_night.mat');
fday_lt = ((fday - 5/24) - floor (fday - 5/24)) * 24;

k = ones([1,length(f107)])';
fday_f107 = datenum(1970.*k', 1.*k', 1.*k', 0.*k', 0.*k', 0.*k'+f107(:,2)');



L = f107(:,3) == 32767;
f107(L,3) = 0.0;
% moving average : for EUVAC
windowSize = 81;
F = ones(1,windowSize)/windowSize;
euvac = conv(f107(:,3),F);
%# if you didn't use 'valid', Df is larger than D. To correct:
halfSize = floor(windowSize/2);
euvac = euvac(halfSize+1:end-halfSize);


fid = fopen('/data/backup/mnair/longp/jic_isr_fejer.txt','wt');

for i = 1:length(fday),
    
    k = find(floor(fday_f107) == floor(fday(i)) );
    f107d = euvac(k);
    [y,m,d] = datevec(floor(fday(i)));
    doy = dayofyear(y,m,d);
    
    fprintf(fid,'%3d %6.3f %6.3f \n', doy, fday_lt(i), f107d);
end;


%% load the fejer output\
% The processing was done on the wondows computer
fejer= load('/data/backup/mnair/longp/jic_isr_fejer.out');

% now creating a data matrix similar to that of eef supplied by
% alken. Note that some columns with 0s will be there
% mutiplying the drift with 23000 nT for getting eef

eef(:,1) = fday - datenum(2000,1,1);
eef(:,4) = fday_lt;
eef(:,6) = drift*23000/1e6;
eef(:,15) = fejer(:,2)*23000/1e6;


%% 
%% Remove the climatological variation from CHAMP EEF data
% writing a file as input to Patrick's eefm model


% load f107 data 
load /data/backup/mnair/indices/f107.dat;
load /data/backup/mnair/longp/eef_data_2000_2009.mat eef;
% # Field 1: timestamp (UT)
% # Field 2: longitude (degrees)
% # Field 3: latitude (degrees)
% # Field 4: local time (hours)
% # Field 5: season (day of year)
% # Field 6: eastward electric field (mV/m)
% # Field 7: equatorial vertical electric field at 105 km altitude (mV/m)
% # Field 8: equatorial (UxB)_z 105 km altitude (mV/m)
% # Field 9: DC current shift (A/m)
% # Field 10: CHAMP peak current value (A/m)
% # Field 11: KP
% # Field 12: F10.7 (W/m^2)
% # Field 13: F10.7A (W/m^2)
% # Field 14: R^2 (coefficient of determination)
% # Field 15: R (correlation of CHAMP and model profiles)
% # Field 16: chi^2
% 

k = ones([1,length(f107)])';
fday_f107 = datenum(1970.*k', 1.*k', 1.*k', 0.*k', 0.*k', 0.*k'+f107(:,2)');

L = f107(:,3) == 32767;
f107(L,3) = 0.0;
% moving average : for EUVAC
% EUVAC = 0.5*(F10.7+F10.7A), where F10.7A is the 81-day moving average of F10.7
windowSize = 81;
F = ones(1,windowSize)/windowSize;
dummy = conv(f107(:,3),F);
%# if you didn't use 'valid', Df is larger than D. To correct:
halfSize = floor(windowSize/2);
f107a =dummy(halfSize+1:end-halfSize);
euvac = 0.5 * (f107a + f107(:,3));

unix_time_stamp = int32((eef(:,1) - datenum(1970,1,1))*24*3600);

fid = fopen('/data/backup/mnair/longp/eefm_model_input.txt','wt');

for i = 1:length(eef),
    
    k = find(floor(fday_f107) == floor(eef(i,1)) );
    f107d = euvac(k);
    [y,m,d] = datevec(floor(eef(i,1)));
    doy = dayofyear(y,m,d);
    
    fprintf(fid,'%8.3f %3d %6.4f %d %6.3f \n', eef(i,2), doy, eef(i,4), unix_time_stamp(i), f107d);
    euvaceef(i) = f107d;
end;
fclose all;


% to find out if there is any LT dependency for sigma error etc

nbin = 24; % 24 LT hour bins
sigma_ut = zeros([1,nbin]);
data_n = zeros([1,nbin]);
for i = 1:length(ntimestamps),
    bin_id = floor((ntimestamps(i) - floor(ntimestamps(i)))*nbin) + 1;
    if ~isnan(jicamarca_sigma_gt_500(i)),
    sigma_ut(bin_id) = sigma_ut(bin_id) + jicamarca_sigma_gt_500(i);
    data_n(bin_id) = data_n(bin_id) + 1;
    end;
end;

    

%%
L =fday_jic >= datenum(2001,12,12) & fday_jic <= datenum(2001,12,13);
un_ts = fday_jic(L);
kk = unique(un_ts);
clear drift_all
clear drift_ts
 for i = 1:length(kk),
L = fday_jic == kk(i);
driftdata = jicamarca_isr(L,7);
sigma = jicamarca_isr(L,8);
LL =  sigma < 2;

drift_ts(i) = nanmean(driftdata(LL));
LL = jicamarca_isr(L,6) >=0;
drift_all(i) =  nanmean(driftdata(LL));
end;
